/* Polynomials.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifndef _polynomials_
#define _polynomials_
#include <boost/operators.hpp>
#include <src/miscellaneous\io\IOFacilities.h>
#include <src/data/Operators.h>

/*! \brief Polynomials class with full algebra set of operations.
	Polynomials class with full of operations : (+,-,*,<,<< (shift operator) )
*/
template<class data = OS_double, class container = particle<OS_double> > 
	struct polynomial : private      
		boost::addable< polynomial <data,container>          // polynomial + polynomial
    ,	boost::subtractable< polynomial <data,container>     // polynomial - polynomial
	,	boost::multipliable < polynomial <data,container>
	,	boost::multipliable2 < polynomial <data,container>, typename container::value_type
	,   boost::less_than_comparable < polynomial <data,container>
	,	boost::shiftable2<polynomial <data,container>, size_t
	> > > > > >,
	public std::unary_function<const data&, data>,
	public derivable< data>,
	public integrable< data >
{
public:
	typedef polynomial									type;
	typedef typename	container::iterator				iterator;
	typedef typename	container::const_iterator		const_iterator;
	typedef typename	container::value_type			value_type;
	polynomial() {};
	polynomial(const polynomial & p) : value_(p.value_){};
	template <class iterator> 
	polynomial(const iterator & begin, const iterator & end) {
		value_.resize(end-begin);
		std::copy(begin,end,value_.begin());
	};
	inline void resize(const polynomial &from) {
	   if (size() < from.size()) {
		   container temp(value_);
		   value_.resize(from.size());
		   for (OS_size i = 0; i < size(); ++i) {
		    if (i < temp.size())
				value_[i] = temp[i];
			else
				value_[i] = 0;
		   };
	   };
	}
    inline polynomial& operator+=(const polynomial &from)  {
		resize(from);
		iterator Ptr_ = value_.begin();
		const_iterator Ptr_From = from.value_.begin();
		iterator End_ = value_.end();
		const_iterator End_From = from.value_.end();
		for (; Ptr_ < End_ && Ptr_From < End_From; ++Ptr_, ++Ptr_From)
				(*Ptr_) += (*Ptr_From);
		return *this;
    };
   inline polynomial& operator-=(const polynomial &from)  {
		resize(from);
		iterator Ptr_ = value_.begin();
		const_iterator Ptr_From = from.value_.begin();
		iterator End_ = value_.end();
		const_iterator End_From = from.value_.end();
		for (; Ptr_ < End_ && Ptr_From < End_From; ++Ptr_, ++Ptr_From)
				(*Ptr_) -= (*Ptr_From);
		return *this;
    };

	inline polynomial& operator = (const polynomial & rhs) {
		value_ = rhs.value_;
		return *this;
	};
	inline polynomial& operator*=(const polynomial & rhs)  {
		polynomial temp(rhs);
		polynomial my_temp( *this );
		value_.resize(degree() + rhs.degree() +1);
		value_ = 0;
		while (temp.size() > 0) {
			if (temp[0]) {
				(*this) += my_temp*temp[0]; 
			}
			temp >>=1;
			my_temp <<=1;
		}
		return *this;
	};
	polynomial(const value_type &c) {
		value_.resize(1);
		value_[0] = c;
	};
	inline polynomial& operator*=(const value_type & rhs)  {
		value_ *= rhs;
		return *this;
	};
	virtual data operator()(const data & x) const { 
		
/*
we can not use Horner method for evaluation due to boost interval algebra mismatch (for an intervall I, I*I is not equal to pow(I,2) ).
		const data & a = x;
		OS_size degree_ = degree(); 
		data y = ( const data& ) value_[degree_];
		OS_size i = degree_;
		do {
			--i;
			y = y * a + ( const data& ) value_[i];
		} while (i != 0);
		return y;
		
we use pow instead (thus this is a slower evaluation) : 
*/

		const data & a = x;
		OS_size degree_ = degree(); 
		data y = ( const data& ) value_[0];
		for (OS_size i = 1; i<= degree_; ++i) {
			y += value_[i]*pow(a, (OS_int) i);
		}
		return y;

	};

	inline OS_size degree() const {
		OS_size deg;
		for (deg = value_.size()-1; deg >=0; --deg) {
			if (value_[deg] != 0) break;
		};
		return deg;
	};
	inline polynomial& operator>>=(OS_size _Pos) {
			polynomial from(*this);
			value_.resize(size() - _Pos);
			value_ = 0.;
			iterator Ptr_ = value_.begin();
			iterator Ptr_From = from.value_.begin()+_Pos;
			iterator End_ = value_.end();
			for (; Ptr_ < End_; ++Ptr_, ++Ptr_From)
					(*Ptr_) = (*Ptr_From);
			return *this;
	};
	inline polynomial& operator<<=(OS_size _Pos) {
			polynomial from(*this);
			value_.resize(size() + _Pos);
			value_ = 0.;
			iterator Ptr_ = value_.begin()+_Pos;
			iterator Ptr_From = from.value_.begin();
			iterator End_ = value_.end();
			for (; Ptr_ < End_; ++Ptr_, ++Ptr_From)
					(*Ptr_) = (*Ptr_From);
			return *this;
	};

	friend inline OS_bool operator < (const polynomial &lhs, const polynomial &rhs) {
		OS_size degl = lhs.degree(); 
		OS_size degr = rhs.degree();
		if (degl != degr )
			return degl < degr;
		for (OS_size deg = degr; deg >=0; --deg) {
			if (lhs.value_[deg] != rhs.value_[deg])
				return (lhs.value_[deg] != rhs.value_[deg]);
		}
		return true; //equality
	}
	friend inline std::ostream& operator<< (std::ostream& o, const polynomial & t){
		o<<t.value_;
		return o;
	};
	inline void operator= (container lfs){
		value_ = lfs;
	};
	inline value_type& operator[](OS_size _Pos) const {	
		return value_[_Pos];
	};
	inline OS_size size() const {	
		return value_.size();
	};
	virtual smart_ptr_function get_derivative() const {
		type* result( new type ( (*this) >> 1) );
		OS_size deg = result->degree();
		for (OS_size i =0; i <= deg ; ++i) {
			( (*result)[i]) *= (i+1);
		};
		return smart_ptr_function(result);
	};
	virtual smart_ptr_function get_integral() const {
		type* result( new type ( (*this) << 1) );
		OS_size deg = result->degree();
		for (OS_size i =1; i <= deg ; ++i) {
			( (*result)[i]) /= i;
		};
		return smart_ptr_function(result);
	};
	container value_;
};


template <class polynome >
struct polynomial_canonical_base {
	typedef typename polynome::value_type value_type;
	polynome operator[](OS_size deg){
		polynome result( (value_type) 1);
		result <<= deg;
		return result;
	}
};

#endif // _polynomials_

